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ABSTRACT 

Recently, the existence of geometrically thick dust stractures in Active Galactic Nuclei (AGN) 
has been directly proven with the help of interferometric methods in the mid-infrared. The ob- 
servations are consistent with a two-component model made up of a geometrically thin and 
warm central disk, surrounded by a colder, fluffy torus component. Within the framework of 
an exploratory study, we investigate one possible physical mechanism, which could produce 
such a structure, namely the effect of stellar feedback from a young nuclear star cluster on 
the interstellar medium in centres of AGN. The model is realised by numerical simulations 
with the help of the hydrodynamics code TRAMP. We follow the evolution of the interstellar 
medium by taking discrete mass loss and energy ejection due to stellar processes, as well as 
optically thin radiative cooling into account. In a post-processing step, we calculate observ- 
able quantities like spectral energy distributions and surface brightness distributions with the 
help of the radiative transfer code MC3D. The interplay between injection of mass, super- 
nova explosions and radiative cooling leads to a two-component structure made up of a cold 
geometrically thin, but optically thick and very turbulent disk residing in the vicinity of the 
angular momentum barrier, surrounded by a filamentary structure. The latter consists of cold 
long radial filaments flowing towards the disk and a hot tenuous medium in between, which 
shows both inwards and outwards directed motions. With the help of this modelling, we are 
able to reproduce the range of observed neutral hydrogen column densities of a sample of 
Seyfert galaxies as well as the relation between them and the strength of the silicate 10 |J.m 
spectral feature. Despite being quite crude, our mean Seyfert galaxy model is even able to 
describe the SEDs of two intermediate type Seyfert galaxies observed with the Spitzer Space 
Telescope. 

Key words: Galaxies: nuclei - Galaxies: Seyfert - ISM: dust, extinction - Radiative transfer 
- Hydrodynamics - ISM: evolution. 



1 INTRODUCTION 

Within the so-ca l led Unified Scheme of Active Galactic Nuclei 
( lAntonucci|[l993l : lUrrv & Padovanilll995h . a torus-like geometri- 
cally thick gas and dust component is coined to explain two ob- 
served classes of Active Galactic Nuclei (AGN) by one intrinsically 
unique AGN model. This torus is illuminated by the accretion disk, 
surrounding the central black hole. As the dust opacity peaks in the 
UV/optical wavelength range, where also the accretion disk emits 
maximally, most of its light is absorbed and reemitted in the form of 
a pronounced peak in the mid-infrared Zanders et al.ll989l) . which 
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is frequently observed in AGN. This gives rise to the separation 
into type 2 sources, where the torus is viewed edge-on and type 1 
sources, where a face-on view allows the visibility of the accretion 
disk, which shows up in form of the so-called Big Blue Bump in the 
spectral energy distribution (SEP ) of AGN. This geometrical unifi- 
cation idea was first proposed bv lAntonucci & Miller. (1985i) , after 
broad emission lines - arising from fast moving gas close to the 
central black hole (the so-called Broad Line Region) - have been 
detected in polarised light in the Seyfert 2 galaxy NGC 1068, which 
was explained by scattering by dust and electrons above the torus 
opening. In direct light, type 2 sources only exhibit narrow emis- 
sion lines, originating from slower gas motions further out, where 
they are less affected by the gravitational potential. This is the so- 
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called Narrow Line Region, located beyond the opening of the torus 
funnel. The first direct evidence for a dusty torus came from in- 
terferometric observations with the help of the MID-infrared In- 
terferometer (MIDI), which revealed geometrically thick dust dis- 
tributions on par s ec scale in the two Seyfert galaxies N GC 1068 
jjaffe et alJIlOOl: IPoncelet et al.ll20q6l : iRaban et al.ll2008l) and the 
Circinus galaxy i Tristram et alfcOOTir 

Up to now, no conclusive, physical model for the distribution 
of gas and dust within these tori exists. This is due to the fact that 
very complicated and yet not fully understood astrophysical pro- 
cesses are thought to happen in the centres of AGN. Additionally, 
their large distances do not allow direct imaging observations. The 
major problem for the persistence of geometrically thick gas and 
dust distributions is to obtain stability of the vertical scale height 
against gravity. So far, se veral models have been put forward: 
iKrolik & Begelmanl l ll988h proposed that the torus is made up of 
clumps which possess supersonic random velocities, maintained by 
transferring orbital shear energy with the help of sufficiently elasti c 
collisions between the clumps (see also lBeckert & Duschlll2004l) . 
To reach this elasticity, high magnetic field strengths are neces- 
sary. As these tori are made up of a multiphase mixture of gas and 
dust with cold, warm and hot components, a clumpy structure also 
helps to prevent the dust from being destroyed by hot surround- 
ing gas. Further evidence for dumpiness or a filamentary structure 
of the cold component - in this case mainly for the distribution 
of neutral gas - comes from X-ray measurement s of the absorbing 
column density distribution ( iRisaliti et ai]|2002h and their combi- 
nation with measurements of the strength of a spectral feature of 
a dust ingredient (se e also discussion in Sect. |6.3I >. Very recently, 
iTristram et al.l ^20071) found hints for dumpiness in the dust torus of 
the Circinus galaxy by means of interf erometric observations in the 
mid-infrared. IWada & NormanI ( l2002h performed hydrodynamical 
simulations of the interstellar medium in centres of active galaxies 
under the assumption of starburst conditions. A very high super- 
nova type II rate in a narrow sheet around the midplane (follow- 
ing in situ star formation in the densest region) is able to puff up 
an initially rotationally supported thin disk. In an alternative sce- 
nario, no torus is needed, but the necessary obscuration is given 
by du sty clouds, which are embedded into a hydromagnetic disk 
wind l lKoni gl & Kartie 1 9941 and discuss i on in ^Elitzur 200l). The 
m ost recent approach c omes from lKrolikI ( |2007|) . following an idea 
of |Pier&Kroii3 ( ll992l) . His idealised analytical calculations show 
that the scale-height of AGN tori can be stabilised against gravity 
with the help of infrared radiation pressure. 

A second approach towards modelling gas and dust struc- 
tures of AGN is via radiative transfer calculations of specific 
dust distributions, which are motivated by si mplified astrophysical 
scenarios. The most up-to-date simulatio ns i Nenkova et alj 200^ 



iHonig et al.l2006l ; ISchartmann et al.l2003 : lNenkova et alj2008j |bl) 
feature clumpy dust distributions. From these calculations, we get 
some rough idea of the geometrical distribution of dust within 
AGN from comparison with high resolution spectral (e.g. NACO 
or Spitzer) as well as interferometric observations with MIDI (the 
MID-infrared Interferometer). 

The goal of this work is to explore the effects of stellar feed- 
back from a nuclear stellar duster on the formation and evolution of 
AGN tori in terms of an exploratory study. In this scenario, the torus 
is mainly built up by gas loss due to stellar evolution. Energetic su- 
pernova explosions provide an effective structuring mechanism. To 
this purpose, we deploy three-dimensional hydrodynamic simula- 
tions. Eventually, this model should account for the obscuration as 
well as the feeding of the central source. In a second step, we con- 



vert the gas density distribution into a dust distribution, which we 
then use for continuum radiative transfer calculations, which can 
later be compared to observational results. In this work, we focus 
on low-luminosity AGN (Seyfert galaxies), as they are more abun- 
dant in the local universe and, therefore, their nuclear regions are 
accessible via current generation interferometric instruments. 

After having discussed the main physical assumptions of our 
work in Sect. (2] we will briefly present the numerical method we 
use (Sect.O. In Sect.|4l we describe the results of this exploratory 
study in terms of the evolution of our standard torus model. Further- 
more, we characterise its final state and explain the results from a 
study of varying several model parameters. After a critical discus- 
sion (Sect. |5}, we compare our torus simulations to recent obser- 
vations of Seyfert galaxies (Sect.|6ll in terms of the spectral energy 
distributions of dust re-emission, as well as gas extinction column 
densities and summarise our work (Sect.[7](. 



2 PHYSICAL ASSUMPTIONS OF OUR MODEL 

The main assumption of our model is that a young star cluster ex- 
ists in those nuclei of AGN, which possess a torus. Evidence for this 
was given recently by observations of a nu mber of Seyfert gal axies 
observed with adaptive optics techniques jDavies et al.ll2007h and 
for the special case of NGC 1068 wit h the help of HST/NICMOS 
images l lGallimore & Matthewsll2003l and references therein), one 
of the closest - and therefore best studied - Seyfert 2 galax- 
ies. We model this nuc lear stellar cluster in form of a Plummer 
profile ( |Plumme3ll91ll) of the gravitational potential and assume 
that it was built up durin g a short-duration st arbursQ in concor- 
dance with the findings of lOavies et alj ( l2007h . which forms a so- 
called Coeval Stellar Population (CSP). The first phase of evolu- 
tion of the stellar population is expected to be very violent, such 
that supernova type II explosions might evacuate most of the nu- 
clear region from gas and dust. During this phase, we do not ex- 
pect that a stable gas and dust distribution (torus) exists. Accord- 
ing to our population synthesis modelling with the st arburst 99- 
code (Leitherer et al. 1999; Va zquez & Leitherejr2005h . this vio- 
lent phase lasts for approximately 40 million years. Within about 
the same time period, double systems of lower mass stars have been 
able to form. According to mod els which include binary systems 
jde Donder & Vanbeverenir2003h . the rate of supernova type la ex- 
plosions - arising from so-called double degenerate systems - be- 
gins to rise steeply and starts to dominate the input of energy. At the 
same time, mass injection should be dominated by the emission of 
planetary ne bulae from st ars in the intermediate mass range (from 
1.5 to 8 Mg. lKwokll200i) . 



2.1 Mass input 

Stars of all masses lose material during their lifetime in form of 
winds of different strengths, mainly happening within short periods 
(e.g. at the tip of the Red Giant Branch). For the case of high-mass 
stars, also supernova type II explosions have to be taken into ac- 
count. The time dependent mass loss rate of a CSP can in principle 
be calculated with the following integration over the initial stellar 
mass for a given (time dependent) mass loss rate of individual stars 
m{t, rrii) and initial mass function (IMF) ipljrii): 



We assume a starburst duration of 40 Myrs in the following. 
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M{t) = 



rh{t, mi) ijjijni) dmi 



(1) 



As both functions are not very well known, a number of sim- 
plifying assumptions have to be applied. We use the model of 
ljungwiert et al.l l lioOll) . which summarises the various phases of 
mass loss in one delta peak for the intermediate mass stars (IMS) 
and the high mass stars (HMS) and two delta peaks for low mass 
stars (LMS). For the IMS, the delta peak represents the AGB 
{Asymptotic Giant Branch) wind peak and for the HMS the SNe 
explosions plus winds. For the LMS, it represents the wind loss 
at the tips of the RGB {Red Giant B ranch) and AGB. The IMF is 
modelled according to lScald l ll998l) as a broken power law. Using 
the stellar initial-final mass relation and the stellar lifetime-mass 
relation f rom the Padua s tellar evolutionary models jBressan et al.l 
M arigo et al. I [1993), an approximate, normalised mass loss 
rate function can be derived: 



Mn{t) 



5.55 • 10" 



f + 5.04- 106 yr' 



(2) 



where the subscript n denotes that it is normalised to the ini- 
tial mass of the CSP. This leads to a mass loss rate of 1.2 • 
10~^ Mq yr~^ ^'^p)^ (normalised to a total stellar mass of 1 A'/q) 
at the beginning of our simulations (after 40 Myr). Recent ob- 
servations tend to find rather top heavy IMFs in the centres of 
galaxies. This means that the number of massive stars is increased 
compared to low-mass stars. It has been shown for the centre of 



the Milky Way jStolte et alj l200l |2005| ; iNayakshin etal 
IPaumard et al.ll2006l) . the Androm eda galaxy JBender et al 
as well as more distant galaxies l lScalolll990| - |Elmegreenl 
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Therefore, we use a five times larger value of Mn = 6.0 • 
10"^ Mq yr-^ in our standard model, which is kept con- 
stant over the whole evolution of the simulation. For simplification 
reasons, we assume that all of the mass is injected in form of plan- 
etary nebulae. 

Numerically, mass injection means a source term in the conti- 
nuity equation as well as the momentum equation, when taking the 
ejection velocity of the gas into account. According to the mass in- 
jection rate, the number of planetary nebulae (PNe) to injec(f| into 
the whole model space is calculated. Residual mass is injected into 
the model space within the subsequent time step. The position of 
the ejecting star is chosen randomly following the distribution of 
stellar mass. The ejection velocity of the PN is chosen such that 
it matches the velocity dispersion (assumed to be independent of 
position) plus the rotational velocity of the stars in the central star 
cluster. As we cannot follow the evolution of single stars explicitly, 
all of the ejected planetary nebulae have the same gas mass mpN, 
which corresponds to the mean mass by weighting the initial-final 
mass relation (given by ^eidemann 200(Q with the Initial Mass 
Function (IMF), which in t his mass regime is well approximated 
by a Salpeter-like function ( ISalDetei|[l95^ . This leads to an av- 
erage mass of 2.2 Mq per planetary nebula, which we distribute 
homogeneously over 3^ cells surrounding the randomly chosen po- 
sition. 

A detailed modelling of the temperature structure of a plan- 
etary nebula is beyond the reach of our current resolution, as we 

As the main mass input at the time of our simulations is given by the 
ejection of PNe, we further refer to this mass input mechanism as the plan- 
etary nebulae injection, although all mass loss effects mentioned above are 
included. 



are not able to resolve the small-scale physical processes. As it will 
finally thermally merge with the surrounding gas on short time- 
scales, we choose the temperature to be the same as the gas in the 
cells where the PNe are injected. 



2.2 Energy input 

Our simulations start after roughly 40 Myr, when the violent phase 
of supernova II explosions (e.g. a 40 Mq mass star lives only 
about 5 Myr) has ended. At roughly the same time, the rate of 
SN la explosions begins to dominate the energy input into the 
ISiy(f| - a result fo und in population synthesis mo delling includ- 
ing binary systems jde Ponder & Vanbeverenll2003h . According to 
today's theoretical understanding, SNIa are thought to be ther- 
monuclear disruptions of accreting white dwarfs, which reach the 
Chandrasekhar mass {^1.4 Mq), being the maximum mass 
which can be supported by electron degeneracy pressure. The de- 
termination of SNIa rates for the considered nuclear stellar sys- 
tem in our simulations is complicated. Observations have to rely 
on low number statistics. Additionally, only the stellar content of 
th e whole galaxy can b e taken into account. A recent publication 
bv lSullivaneraU ( l2006l) describes the SN la rate for young and old 
stellar populations and finds that it can be well represented by the 
sum of 5.3 ± 1.1 ■ 10"^* SNeyr"^MQ^ (old population of stars) 
and3.9±0.7-10"'' SNeyr"^(MQ t/r"^)"^ (young population of 
stars). The rate for the old stellar population is normalised to 1 AIq 
of stars and the latter rate to a star formation rate of 1 Mq yT~^. 

For an assumed duration of the starburst of 40 Myr, this re- 
sults in a normalised SN la rate of roug hly 10~"SNeM-Vr"^ 
for our modelling. Due to the large uncertainties in theoretical 
derivation of supernova rates and as it is doubtful whether the ob- 
servational derivations are appropriate for our simulations, we con- 
sider the SN la rate a more or less free parameter within our mod- 
els. In_an_alternatiyestarb^ scenario for nuclei of Seyfert galax- 
ies, IWada & Norma j ( |2002|) use a much larger rate of roughly 100 
times our SN-rate, namely I SN detonation per year, restricted to 
a narrow strap around the midplane of their model space. This is 
motivated by violent star formation within a dense gas disk in the 
galactic midplane, accounting for SN II explosions. After a SN-rate 
parameter study, a value of 10^^'' SNe Mq^ yr~^ was chosen for 
our standard model. It is assumed to be constant in our simulations, 
as our current period of time evolution (1.2 Myrs) is short compared 
to the expected evolutionary time scales of AGN tori, of the order 
of 100 Myrs. From the supernova rate for the whole model space, 
the number of supemovae exploding in each timestep can be deter- 
mined. When a supernova is launched, the position is determined 
randomly, but according to the underlying stellar distribution. The 
energy is injected into a single cell in the form of thermal energy. 
Additionally, a mass of 1.4 Mq is added to this cell, also ensuring 
numerical stability. We use an efficiency of r;sN = 0.1, reducing 
the input of thermal energy to 10''° erg. This factor accounts for 
energy losses in the early phases of evolution of the supernovae, 
which we cannot model at the current resolution. The problem of 
radiating away most of the introduced thermal energy due to vigor- 
ous cooling in very dense regions before a supernova shell can form 
is mitigated by introducing a cooling delay of several timesteps. Di- 
rect modelling of an expanding shell as initial condition for a super- 
nova explosion is not feasible at the resolution of our current sim- 



This is only the case, if the so-called double-degenerate scenario of two 
merging white dwarfs is realised in nature. 
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Table 1. Parameters of our standard model. 




10' 10^ 10' 10" 10^ 10*^ 10^ 10" 10^ io'° 
T [K] 



Figure 1. Cooling curve jPlewdl 19951) . The hatched regions denote three 
temperature regimes: cold neutral medium (CNM), warm neutral medium 
(WNM) and hot ionised medium (HIM). The thick red line segments mark 
the thermally stable parts of the cooling curve. Dotted lines refer to cooling 
curves used by other groups, see text. 



ulations. These shock fronts efficiently dissipate their kinetic en- 
ergy when interacting with dense blobs or filaments or dense shock 
fronts of neighbouring supernova explosions, transferring their ki- 
netic energy back to thermal energy. Similar procedures of energy 
input have b een applied to many diff erent astrophysical scenarios. 
For example IWada & NormanI ( |2002|) also introduce the complete 
amount of supernova energy into a single cell in form of thermal 
energy within a starbursting gaseous disk with even higher ambi- 
ent densities compar ed to our simulations. Pur ely thermal energy 
is also introduced byE oung &MacLowll2006l) for the case of the 
interstellar medium in our own galaxy, where the mean gas density 
is much lower compared to our computations. Given the aims of 
this exploratory study, we therefore think that our energy injection 
mechanism is an appropriate treatment of this feedback process. 



2.3 Optically thin cooling 

A mixture of neutral and ionised gas mainly cools via the transfor- 
mation of kinetic energy to radiation by collisional processes. This 
also means that the efficiency of cooling is a sensitive function of 
the constituents of the gas. The radiation can only escape, if the 
medium is optically thin for these wavelengths, otherwise cooling 
is inefficient. 

Fig. [T] sho ws the effecti ve radiative cooling curve applied in 
our simulations jPlewa] 1995h, whic h was calculated with the help 
of the Cloudy code "( lFerlanj[T993h . For the case of our cooling 
curve, a purely collisionally ionised optically thin plasma with solar 
abundances wa s used. 

Following Ijoung & Mac Low! l l200d) , such a gas is thermally 
stable at temperatures, where the slope in logarithmic scale is 
steeper or equal to one. For the case of our cooling function, we 
find five distinct stable regions: 10 - 175 K, 4 800 - 17 400 K, 



1.00 • 10 K. They are marked with the thick red line segments 
in Fig. [T] Overplotted as yellow dotted lines are the often used 
cooling functions from iDalgamo & McCravl ( Il972h - where only 
the part below 10'* K is shown - with ionisation fractions of 0.1 
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Value 


Parameter 


Value 




6.6 • 10'^ Mq 


a 


0.5 


Af* 


1.9 ■ 10^ Mq 


Tini 


2.0 ■ lO'^K 




1.2 • 10* Mq 


Mn 


6.0 ■ lO"'' MQ/{yrMQ) 


Rc 


25 pc 


MpN 


2.2 Mq 




5pc 


SNR 


10-1" SNe/(yr Mq) 


Rout 


50 pc 


r 


5/3 




165 km/s 







Mass of the black hole (A/bh), normalisation constant of the stellar po- 
tential (A/,), initial gas mass (A/g"!,), cluster core radius (-Rc), torus ra- 
dius (Rt), outer radius (i?out), stellar velocity dispersion (cr,), exponent 
of the angular momentum distribution (a), initial gas temperature (Tini), 
normalised mass injection rate (A/n), mass of a single injection (A/pn), 
supernova rate (SNR) and adiabatic exponent (F). 

for the upper curve and 0.01 for the lower one. The violet dot- 
ted line results from equilib rium ionisation cooling, calculated by 
[Sutherland & Dopital l fl993l) for the case of solar metallicity. As can 
be seen from this, the cooling curves are in relatively good agree- 
ment in the high temperature regime (above the ionisation temper- 
ature of hydrogen at around 10* K). At temperatures lower than 
10* K, the values of the assumed cooling curves are most uncer- 
tain, as in this temperature regime, the approximation of thermal 
equilibrium is usually not valid. Therefore, large deviations exist 
between the different calculational methods. It is also a matter of 
debate, which ionisation fraction to assume for the low temperature 
part. We decided to use 0.1 ( black curve and upper yellow curv e 



in Fig. [T} in concordance wi th de Avillez & Breitschwerdj ( |2004|) . 
whereas ljoung & Mac Lowl ( l2006h assume an ionisation fraction 
of 0.01 (lower yellow curve). For temperatures lower than 10 K, 

we set the cooling rate to zero. 

Similar as in I Joung & Mac Low! ( l2006h . we make a division 
into three temperature regimes: the Cold Neutral Medium ( CNM) 
up to 175 K, where the first stable part of the cooling curve ends, 
the Warm Neutral Medium (WNM) in the intermediate temperature 
regime up to the ionisation temperature of hydrogen and the Hot 
Ionised Medium (HIM) for higher temperatures (compare to Fig.[T]l. 



2.4 Parameters of our standard model 

The parameters we use for our standard model are summarised in 
Table[T] In the past subsections, we already motivated the super- 
nova as well as the mass injection rate. For the remaining parame- 
ters, short explanations are given in this paragraph. 

The underlying gravitational potential is made up of two com- 
ponents: 

D The potential of the nuclear star cluster, which is modelled ac- 
cording to a Plummer distribution: 



+ Ki 



(3) 



The Newtonian point-like potential of the nuclear black hole 

GMbh 



6bh(7') — 



(4) 



A core radius Rc=25 pc, a black hole mass A/bh = 6.6 • 
10^ Mq and a normalisation factor for the stellar distribution of 



1.9 • 10 Mq is used. The latter leads to a total integrated 
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stellar mass within the core radius of 6.7 ■ 10 Mq. Our initial 
c ondition comprises of a 7 TM-model similar to the one discussed 
in lSchartmann et 2A\ ( l2005h . but with a Plummer potential (as given 
above) and the parameters listed in Table|T] The angular momentum 
of the stellar content, the gas in the initial condition as well as the 
injected gas is given by the following distribution: 

ispcc = VG(A/bh +A/.(J?T))i?T {-§^Y ■ (5) 

with a constant a = 0.5 and a torus radius of _Rt = 5pc, where 
the torus radius is defined as the minimum of the effective potential. 
The turbulent pressure of the TTM-model is substituted by thermal 
pressure. For the case of our standard model, the velocity dispersion 
of the assumed dust clouds in the TTM-model of cr, — 165 km/s 
then corresponds to an initial temperature of Tini = 2.0 • 10® K. The 
mass loss rate was chosen to be 6.0 • 10"^ MQ/{yr Mq), which 
means a mass injection of 6.0 • 10^'"* solar masses of gas per year, 
normalised to one solar mass of stars (see discussion in Sect. 12. lb . 
A SN-rate of 10^^" supernovae per year and normalised to 1 A/q 
in stellar mass was used (see Sect. |2.2| for further discussion). We 
apply the equation of ideal gas as a closing condition with an adi- 
abatic index of an ideal monoatomic gas (F = |). The gas has so- 
lar abundances in our simulations, with a mean molecular weight, 
which changes from 0.6 (fully ionised gas with solar abundances) 
to 1.2, following a smooth transition at the ionisation temperature 
of hydrogen. Most of the simulations discussed in this paper ran 
for an evolutionary period of 10 orbits at the torus radius {global 
corresponding to 1.2 ■ lO'' yr. Note that this is a tiny frac- 
tion of the age of the stellar population. Therefore, evolutionary 
effects can be ignored. 

3 NUMERICAL METHOD 

For our si mulations, we use the 3D radiation hydrody namics code 
TRAMP l lKlevll 19891 ; fKlahjl 19981 : [klahretalJ[T99^. which uses 
similar numerical scheme s as the Zeus code l lStone & NormanI 
ll992iJ lbl: [stone et al.iri992h . It solves the equations of ideal hydro- 
dynamics with the help of an Operator Splitting technique on a 
spherical grid that is logarithmically spaced in the radial direction. 
Using this method, one has always an optimal resolution and at 
the same time more or less boxy e. g. equi-spaced grid cells. For 
the ad vection step, the monotonic transport scheme of Ivan Lee3 
( 1 19771) is applied, which is a spatially second-order-accurate up- 
wind method. The radiation transport in TRAMP is treated in 
the fl ux-limited diffusion approximation using the flux limiters of 
iLeve rmore & Poimaning (1981). Additional source terms for mass 
input (Sect. 12. lb and energy input (Sect. 12. 2b are treated in a sep- 
arate step. TRAMP discretises the hydrodynamic equations in the 
finite volume description and, therefore, locally and globally con- 
serves advected quantities. The force step is calculated using a 
finite difference discretisation. Shocks are treated by introducing 
a non-linear artificial viscosity, which smears out discontinuities 
over several cells and, thereby, results in the correct generation of 
entropy in the shock, the correct shock velocities and suppresses 
post-shock oscillations. It is treated as an additional contribution to 
the pressure. The non-linearity ensures that the artificial viscosity 
is large in shocks, but negligible elsewhere. 



3.1 Domain decomposition and boundary conditions 

In order to gain resolution, we split the whole computational do- 
main into three parts: the inner one ranges from 0.2 to 2.0 pc, 
the middle one from 1.0 to 10.0 pc and the outer part from 5.0 
to 50.0 pc (logarithmic grid). Towards the centre, the size of the 
single cells decreases and, additionally, the obtained velocities in- 
crease due to the presence of the central gravitational potential and 
the conservation of angular momentum of infalling material. This 
leads to ever smaller timesteps, in order to prevent information to 
be transported along more than one cell during a single timestep 
{Courant-Friedrichs-Levy condition, ICourant & Friedrichsl Il94a : 
iRitchmver & Mortonll967h . Hence, we are able to calculate a much 
longer time period for the outer domain with the same amount of 
CPU usage and the same grid size, compared to the inner domain. 
On the other hand, further in, processes act on much shorter time- 
scales and only a shorter time evolution is needed in order to ob- 
tain a steady state. Therefore, we calculate 10 global orbits for the 
outer domain, but only the last orbit for the middle domain and one 
tenth of an orbit for the innermost part. This leads to a comparable 
amount of CPU time of approximately one weak per domain. 

In the first step, the outer domain is calculated. The physical 
state of material at the position of the outer boundary of the mid- 
dle model is stored in a file every 0. 1 % of a global orbit. Outflow 
boundary conditions are used radially inwards as well as outwards. 
This means we allow for the outflow of matter from the grid and 
minimize the reflection of waves at the grid edge, while not allow- 
ing for inflow, thus preventing the boundary cells from generating 
a numerical instability. During the calculation of the middle do- 
main in the second step, the inner radial boundary is set to outflow, 
whereas the outer boundary is set by the afore stored quantitiejfl 
which can result in in- or outflowing material, thanks to the wide 
overlap of the domains. At the same time, physical quantities ad- 
vected across the position of the outer boundary of the inner domain 
are stored. In the last step, the same is done for the innermost model 
space. For the final quantitative analysis and the subsequent radia- 
tive transfer calculations, the three domains are combined again. 
This procedure is a reasonable approximation, as in most of the 
simulations we observe infall of gas towards the inner region. Out- 
flowing material can partly be taken into account, as the domains 
have a quite large overlapping region. 

Concerning the other coordinates, we model a range of 90° 
with periodic boundary conditions in 4> direction. In ^-direction, 
where we model from 4 degrees to 176 degrees, as well as for the 
inner- and outermost boundary in radial direction, outflow bound- 
ary conditions are applied in a sense that mass is allowed to flow 
out of the integration domain, but no material can get back into it. 
This means that the values of the ghost cells are set such, that 
the normal derivatives of the boundary values equal to zero for all 
quantities, presenting a smooth continuation of the flow through 
the boundary. Furthermore, we only allow outward directed mo- 
tion, otherwise the normal velocity components are set to zero. For 
each of the domains, we are able to calculate 75 cells in radial, 94 
cells in 6 and 49 cells in 4> direction, leading to cell sizes between 
0.006 pc in the innermost part of the whole domain and 1.5 pc at 
the outer edge. 



* Global orbit always refers to a Keplerian orbit, measured at the torus ^ A linear interpolation in time is used. 

radius at /?x = 5 pc and corresponds to a real time of approximately 1.2 ■ ^ Ghost cells are additional boundary cells, not belonging to the simulated 

10'' yr physical domain, but necessary to implement boundary conditions. 
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3.2 Initial condition 

A hydrodynamically stable initial condition for the case of our 
effective potential - made up of the nuclear stellar cluster, the 
black hole and the additional centrifugal potential due to or- 
bital r otation - is given by the isothermal Turbul ent Torus Model 
(TTM. ICamenzindlll995l ISchartmann et al.ll2005l) . A conical vol- 
ume around the rotation axis with a half opening angle of 4° needs 
to be cut out, because the funnel region of this field configuration 
is dynamically unstable and would cause unphysical behaviour of 
the state vector. 



4 RESULTS 

4.1 Tiie evolution of our standard model 

Within the first few percent of a global orbit, single injected blobs 
of gas are visible, mainly within the central 25 pc, the core radius 
of the stellar distribution (Fig. [2^). Spherical cavities are formed 
due to the explosion of supemovae. While these explosions pre- 
dominantly form shock fronts, the planetary nebulae interact in 
two different w ays, as we could se e in separate simulations of sin- 
gle events (see ISchartmann|[2007h : In the subsonic case (e.g. hot 
surrounding medium), colliding streams tend to stick together and 
form long filaments in perpendicular direction from the initial di- 
rection. The supersonic interaction of two planetary nebulae (e.g. in 
cold medium) will cause the formation of a common shock front, 
which surrounds the point of interaction. With ongoing evolution, 
the single blobs collide to create a filamentary, net-like structure, 
enhanced and further shaped by supernova explosions (Fig.|2j3 and 
c). This process is also amplified by the onset of cooling instabil- 
ity in the dense regions of the filaments. Cooling also leads to loss 
of thermal pressure and, therefore, material tends to fall towards 
the minimum of the potential well, located at the torus radius (5 pc) 
within the equatorial plane of the torus. Around this radius, the fila- 
mentary structure opens out into a geometrically thin, but optically 
thick and very turbulent disk (Fig.|2ji). 

Thi s process is comparabl e to the so-called galactic fountain 
process jShapiro & Fieldri976r) visible in galactic disks. The more 
tenuous the gas, the more effective are the supernova explosions in 
pushing away material. This mechanism finally leads to the large- 
scale, filamentary structure visible in Fig.[2jl. 

As the implemented optically thin cooling scales proportional 
to the gas density squared, only the dense filaments cool on short 
time-scales, whereas the interclump-medium is effectively heated 
by the energy injection of the supernova explosions. Accordingly, 
the gas velocity can appreciably differ in various temperature 
regimes, as will be discussed in Sect. 14. 31 

In order to ensure that the memory of the initial condition is 
completely lost, we start calculating the middle part of the compu- 
tational domain after an evolution of nine global orbits of the outer 
domain. Already after a fraction of the simulated one global orbit, a 
geometrically relatively thin, but optically very dense disk forms in 
the equatorial plane in the vicinity of the minimum of the potential 
(see Fig.[6]and discussion in Sect. l4.4l >. 

The final state of our standard model after ten global orbits 
is depicted in terms of density, temperature and pressure within 
a meridional plane in Fig. [3] A multiphase medium has formed 
(Sect. 14.21 ). The temperature and the density distribution are com- 
plementary in the sense that high density regions possess low tem- 
perature due to rapid cooling and hot regions possess low density 
as material has been blown away by supernova explosions. It is im- 



portant to note that the pressure distribution (right panel in Fig.O 
is far away from pressure equilibrium, as it spans more than five 
orders of magnitude. 

4.2 Multiptiase medium 

Fig. m shows phase diagrams for the total model space of our stan- 
dard model (all three spatial subregimes) after ten global orbits. 
The hatched regions indicate thermally stable parts of the cooling 
curve (see Sect.|2}. Capital letters (A-F) denote the six bumps in the 
probability density function for gas temperature (upper left panel). 
C and D as well as the small excess E can be directly related to 
the stable regions of the cooling curve, because time spent in each 
temp erature range is pro portional to the inverse of the cooling func- 
tion l lGerola et alJl 19741) . At these temperatures, gas cools compar- 
atively slow and therefore assembles in these regimes. The expla- 
nation for the maximum B is identical. The minimum in-between 
A and B evolves, because gas at these temperatures is mostly lo- 
cated in dense clumps. As cooling depends on the square of the 
density, it can cool quite fast, even if the cooling curve indicates 
thermally stable behaviour. Peak F can be explained by constant 
energy input from supernova explosions, which allows gas to re- 
side even in the thermally unstab le regime (see also discussion in 
Ide Avillez & Breitschwerdll2005l) . Looking at the mass distribution 
of the various temperature phases (lower left panel), a similar curve 
is obtained, but tilted, as most of the mass is concentrated in the 
cold temperature regime (dense filaments and disk), while the hot 
gas fills most of the volume (supernova blown cavities). 

Plotted with respect to the gas density (second column of 
Fig.|4j, three local maxima are visible (G,H,I). Feature I is caused 
by the very dense disk component, as it is dominant in mass, but 
unimportant in volume. This material had enough time to cool to 
temperatures below 175 K and, therefore, conesponds to maximum 
A and B in the left panel. G can be accounted to regions, where su- 
pernova explosions happened recently, contributing to a large vol- 
ume fraction, but only little to the total mass. In the left column, this 
corresponds to feature F. The bump (H) in between arises from ma- 
terial in the large scale filaments, which is in the process of settling 
down and thereby cooling. 

The distribution of mass and volume with respect to pressure 
is shown in the right column of Fig.|4] As already evident from the 
two-dimensional cuts in Fig. [3] no pressure equilibrium is reached. 
This fact is in qualitat ive agreement with recent observations of 
Ijenkins & Tripij ( |2007|) in our local neighbourhood, who find vari- 
ations in the pressure of the diffuse, cold neutral medium in a range 
of 10^ Kcm"^ < P/ks < 10* Kcm"^ by using Cl absorption 
lines, measured in the ultraviolet spectra of roughly 100 hot stars, 
taken from the HST archive. Instead of a single and sharp peak 
(expected for pressure equilibrium), we obtain three local max- 
ima (right panel in Fig. |4ll. The system tries to reach an equilib- 
rium value (located around peak K). Peak L can be assigned to 
regions with recent supernova action, leading to the highest pres- 
sure values of the simulation and filling a large volume. This is 
the same component causing the peaks F and G. We interpret de- 
viations from a global equilibrium pressure value towards small 
pressure (peak J) as a sign of ongoing thermal instability within 
the fast cooling, densest regions of the flow (see also discussion in 
IJoung & Mac Lowll2006l) . This interpretation is backed by Fig. [3] 
where cold and dense regions correspond to low pressure states and 
hot and tenuous regions to high pressure values. 

In terms of volume and mass filling factors, the following pic- 
ture evolves: Most important in mass is the CNM with 60.4% of 
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Figure 2. Temporal evolution of the density in our standard model, given in a meridional slice of the torus. Shown are four stages of the torus at 0.25, 1.0, 5.0 
and roughly 10.0 global orbits. The scaling is logarithmic. 
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Figure 3. State of the torus after roughly ten global orbits. Shown are from left to right the density distribution, the temperature distribution and the pressure 
distribution of a meridional slice thi'ough the 3D data cube. All images are displayed in logarithmic scaling. 



the total mass, mostly residing in the disk component and the inner 
parts of the long filament, which had enough time to cool to tem- 
peratures of the CNM and occupies the smallest fraction in volume 
(only 5.4%). Second most important (22.6% in mass) is the WNM, 
made up of material in the process of cooling down. Only 16.9% of 
the mass is in the HIM, which occupies most of the space (83.3% 
of the volume). Only 11. 3% of the whole model space is occupied 
by the WNM. 

4.3 Dynamical state of the torus 

The build-up of the disk and torus component can also be quan- 
tified by means of radial velocity distributions, as given in Fig. [5] 



The velocity is averaged over all angles 6 and (j) and separated into 
the three temperature components. For an integration time between 
eight and ten global orbits, the data are plotted every 0.05 global 
orbits (single dots). The solid lines are temporal averages for the 
whole range of displayed data. For comparison, we show the free 
fall velocity of a test particle within the total gravitational poten- 
tial of the black hole and the stellar distribution, scaled down by 
a factor of three. It is a remarkable feature of our simulations that 
the slope of the velocity distributions of all three components are 
similar to free fall. 

Clearly, the hot tenuous medium is able to provide the largest 
thermal pressure support, which slows down the infall towards the 
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Figure 4. Phase diagrams for the total model space of our standard model after an evolution of ten global orbits. The hatched regions indicate thermally stable 
regimes according to the cooling curve (see thick red line segments in Fig.[T). 
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Figure 6. Formation of a turbulent, geometrically thin disk close to the midplane. Shown is the logarithm of the density in units of g/cm'^ as a cut through the 
midplane (upper row) and through the central meridional plane (lower row). Time is given in global orbits. Labels are given in pc. 



mass centre. In contrast to this, cold material sinks much faster to 
the bottom of the potential well, within the long, dense filaments. 



At the border line between the phases with different flow 
structures, Kelvin-Helmholtz and Rayleigh-Taylor instabilities are 
expected to occur, but are not visible in our simulations due to the 
lack of sufficient resolution. They are supposed to increase turbu- 
lence and promote further mixing. Due to their high density, the 
filaments are very stable and single supernova explosions are not 
able to disrupt them. In them, the CNM reaches radial velocities of 
almost 30% of the free fall velocity, the WNM approximately 15% 
and the HIM 10%. This is a further indication that energy in form 
of supernovae or due to the planetary nebulae injection is mainly 
dispensed into the hot medium, which is also visible in the larger 
scatter of the radial velocities of the hot component. 



4.4 Formation of a turbulent disk 



As already discussed in Sect. |4.1| in the vicinity of the torus radius, 
a geometrically rather thin, but optically very dense disk forms. 
Fig. |6] shows the temporal evolution of density for four snap-shots 
of the middle domain. Its dynamical evolution is mainly determined 
by mass inflow through the outer radial boundary. This material 
reaches the middle part both in the form of broad streams as well as 
in blobs and filaments of material and gets structured by supernova 
explosions and planetary nebulae injection. As more and more ma- 
terial falls towards the minimum of the potential, a turbulent, fluffy 
disk emerges between the inner boundary and roughly six or seven 
parsec. We note the remarkably filamentary structure of the disk, 
showing strings of material, elongated in azimuthal direction (for 
reasons identical to those already discussed for the outer domain). 
The difference is that the average flow is here directed in azimuthal 
direction and not in radial direction. 



Stellar feedback in AGN tori 9 




Figure 5. Radial velocities of the outer model space between eight and ten 
global orbits, separated into temperature phases. Negative velocities refer 
to infalling motions. The blue line represents the free fall velocity (fff), 
scaled with a factor of 1/3, the dashed line symbolises the location of the 
core radius of the stellar distribution. 



Altogether, the evolution can be summarised as a collapse 
from a geometrically thick distribution (the torus component) to a 
geometrically thin disk. The initial torus is mainly (thermal) pres- 
sure supported and, therefore, collapses due to radiative cooling 
processes into dense filaments, as cooling and heating of the gas is 
locally unbalanced. This finally leads to the formation of the thin 
disk in the vicinity of the minimum of the gravitational potential, 
close to the angular momentum barrier. Thermal pressure does not 
play a role anymore, as the disk has cooled to temperatures of less 
than 100 K and reaches in some regions even our minimal tem- 
perature of 10 K. The scale height is now determined by random 
motions of the filamentary disk and rotation. This hydrodynami- 
cal turbulence leads to angular momentum transfer outwards and 
accretion towards the centre. 



4.5 Parameter studies 

In order to investigate the dependencies of our model, we con- 
duct ed several parame ter studies, which are briefly presented here 
(see lSchartmannllioOTi for a more detailed discussion). 

D Mass and energy input rate: Changing the mass injection rate and 
the supernova rate, we found that a complementary behaviour is 
obtained. Both, with an increasing supernova rate or a decreasing 
mass-loss rate, the inward bound flow is more and more turned into 
an outflow, as demonstrated for example in Fig. |7] This has also 
been shown in ID hydrodynamic al simulations of the interstellar 
medium in elliptical galaxies by iGaibler et all l l2005h . The tran- 
sition radius between inflowing and outflowing hot gas increases 
with increasing mass loss rate or decreasing supernova rate. Inside 
the transition radius, the curves flatten for models with larger mass 
input rate. This is caused by the additional turbulent pressure com- 
ponent, provided by the PNe injection of gas blobs with random 
velocities. 

Cooling rate: Cooling is a sensitive function of metallicity. The 
larger the metallicity, the larger the cooling rate. We varied the 
cooling rate to account for the unknown metallicity of the gas in 
AGN tori. The most striking difference appears in the broadening 
of the filaments, when increasing the cooling rate, as gas of even 
lower densities is able to effectively cool to CNM temperatures. 
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Figure 7. Comparison of radial velocities within the mass loss rate study for 
the HIM between global orbit eight and ten. Lines represent time-averaged 
(global orbits eight to ten) values. 



?) Mass of the central stellar cluster: Varying the (generally un- 
known) mass content in stars is of interest, as it first determines 
the gravitational potential in the outer part beyond the influence of 
the black hole. Second, it changes the amount of supernova input 
and mass input in equal measure. A higher stellar mass produces 
a deeper potential well and therefore accelerates the gas to higher 
velocities and the transition radius between inflow and outflow of 
the hot temperature phase increases. For the case of a higher total 
stellar mass, the evolving structures are more compact and a more 
massive disk is able to form at the angular momentum barrier. 

D Various initial conditions: In order to check the dependence of 
our results on the initial condition, we started one simulation with 
an initially empty model space (possessing a homogeneous density 
distribution with the lower density threshold of our simulations - 
lO"^*" g/cm"^). As expected, this has no visible consequence on the 
evolution of the gaseous structures and a steady state in terms of 
total gas mass is already reached after less than 3 global orbits. 
Also starting with cold (T = 1000 K) initial conditions, where 
the initial torus is supported by turbulent pressure on small scales 
yields identical results. This small-scale turbulence decays within 
a fraction of an orbit and is dissipated into heat. Additionally, the 
supernova input quickly heats the tenuous medium outside-in to- 
wards the minimum temperature at the torus radius. Therefore, af- 
ter a short time, we again get a very similar behaviour as in our 
standard model. 



5 DISCUSSION 

It needs to be emphasized that within our model, the large-scale 
torus component is not a long-living steady state, but it can only 
exist as long as mass and energy is injected from the stellar cluster. 
After switching off energy and mass input, most of the large scale 
structure will disappear within approximately 1.5 global orbits. The 
turbulent disk in the vicinity of the angular momentum barrier is the 
only remaining structure. In reality, this shutdown will happen nat- 
urally due to the evolution of the nuclear star cluster, but on much 
longer time-scales. It means that - within our modelling - the ex- 
istence of tori is closely linked to the existence and evolution of a 
stellar cl uster in the centre of a galaxy. This fits well to the observa- 
tions bv lDavies et al.l ( l2007h . which reveal a link between nuclear 
starbursts and AGN activity with a certain time delay in-between. 



10 M. Schartmann et al. 



Furthermore, it needs to be stressed that the torus scale height in 
our modelhng is mainly determined by the spatial distribution of 
the stars. 

The total simulated time of our runs is of the order of 10 global 
orbits, which corresponds to 1.2 Myr. This is sufficiently long to 
avoid memory effects from the initial condition and at that point, 
a global dynamical equilibrium concerning the total mass and to- 
tal energy is reached for the outer domain. In the middle and in- 
ner domain, the calculated time evolution of one and 0.1 global 
orbits is not long enough to reach a global dynamical equilibrium 
state. However, the technical limitation of this study did not allow 
a longer evolution. 



5.1 Limitations of our simulations and future work 

This paper presents an exploratory study of a self-consistent hydro- 
dynamical model of gaseous and dusty AGN tori. Due to the sim- 
plifications required by the efficiency of our code, some important 
problems and astrophysical effects remain open: 

® Longer time evolution: With our serial code we are limited in 
temporal evolution to about 1.2Myrs, while typical evolutionary 
times of such systems are of the order of 10 to lOOMyrs. There- 
fore, we currently po rt our routines to the parallel PLUTO code 
jMignone et al.ll2007h . 

© Radiation contribution from central source and stars: The effect 
of the radiation of the central accretion disk is twofold: First of 
all, it contributes to the heating of the dust and gas. Second, its ra- 
diation pressure plays an important role in the dynamics of AGN 
dust tori, as was recently shown by,Krolik ( 2007). With the help of 
our detailed radiative transfer calculations, we compared the accel- 
erations due to radiation pressure of the dust reemission with the 
gravitational accelerations in vertical direction for each cell within 
the computational domain of the snapshot after roughly 10 global 
orbits, using a flux mean dust opacity in combination with a grey 
approximation. In this estimation, infrared radiation pressure sig- 
nificantly contributes only close to the midplane in the inner few 
parsecs from the heatin g source, whi ch is in accordance with the 
detailed calculations of Krolild ( l2007h . Within the large scale fila- 
ments, the ratio between the accelerations due to the infrared radia- 
tion pressure and the gravitational acceleration in vertical direction 
amounts to a few percent only, except for the direct illuminated 
parts of the filaments, where it can also exceed the gravitational 
accelerations and will therefore drive an outflow along the torus 
axis, widening the (so far too narrow) opening angle in our simu- 
lations. However, detailed radiative transfer calculations during the 
whole course of the hydrodynamical sim ulations are nece s sary t o 
finally clarify this question. Opposed to Thompson et al.l ( |2005|) . 
who study radiation pressure-supported starburst disks, we are cur- 
rently interested in the phase following an efficient starburst. By 
then, the stellar luminosity has already dropped by a large factor. 
Giving the low stellar luminosity densities derived from the ob- 
servat ions of a Seyfert galaxy sample discussed in lOavies et al.l 
( l2007h , the stellar luminosity seems to be only a small fraction 
of the luminosity of the central accretion disk at this evolution- 
ary stage. Therefore, we currently neglect its contribution to the 
dynamics, as well as to the heating of the dust within our postpro- 
cessing radiative transfer calculations. 

® Magnetic fields: Weak magnetic fields in combination with shear 
due to differential rota tion cause magneto-rotational instability 
tealbus & Hawlevl[l99ll) , which provides a means to transport an- 
gular momentum and thus enables accretion towards the core (that 



is through the inner boundary of our computational domain). Prob- 
ably coexisting strong magnetic fields in a spatially distinct region 
cause hydromagnetic outflows, which might be preferentially di- 
rected along the torus axis and thereby widening the opening angle 
of the torus. 

® Self gravity and star formation: At the current state of our mod- 
elling, self-gravity of the gas is not taken into account. However, 
examining the Jeans-criterion for each cell of our model space, we 
find that part of the inner turbulent disk, as well as a small fraction 
of the filamentary structure is Jeans unstable and will form stars. 
This is consist e nt wit h the findings of iNavakshin et all ( l2006h and 
IPaumard et al.l ( |2006|) . In the snapshot after 10 orbits, the unsta- 
ble mass amounts to 0.3% in volume or 15% in mass respectively. 
However, it should be emphasized again that our purely hydro- 
dynamical code is not capable to treat the accretion flow towards 
the core correctly, leading to a pile-up of cold material in the disk 
which does not describe the steady state density. Subsequently, su- 
pernova type II explosions might inflate or even partly disrupt the 
evolving dense and geometrically thin turbulent disk. However, one 
should note again that detailed simulations including self gravity 
will be subject of our future work. 

® Further mass input: There will also be material transported into 
the simulated domain of the active galactic nucleus from further out 
of the galaxy. Prominent examples are dusty mini-spirals, of which 
some have been discovered by NACO jPrieto et al.ll2005h or HST- 
ACS observations (Hubble Space Telescope - Advanced Camera 
for Surveys). Simulations l iMacieiewskill200^ show that with the 
help of such nuclear spirals, enough gas inflow can be produced in 
order to power luminous nearby AGN. 

© Detailed data comparison: After including part of the discussed 
additional features, one should do a detailed comparison with ob- 
servational data. Starting from the simulated surface brightness dis- 
tributions, we are able to calculate visibilities in order to compare 
to recent mid-infrared interf erometric observations. Furthermore, 
SINFONI observations (e. g. lOavies e t al. 2007) of Seyfert nuclei 
and starbursting galaxies will give us important constraints for our 
modelling, e. g. in terms of the nuclear stellar cluster, as well as 
fuelling of the central active region. 



6 COMPARISON WITH OBSERVATIONS 

The only possibility to directly resolve the dust structure in nu- 
clei of Seyfert galaxies so far is with the help of the MID-infrared 
Interferometric instrument MIDI at the VLTI. T hereby, tori were 
found in the clos est Seyfert nu clei NGC 1068 lljaffe et alj|2004 
iRaban et al.l2008l) and Circinus jTristram et alj2007[) . Remarkably, 
both objects are well modelled by a warm, more or less spheri- 
cal component and a less extended, hotter elongated component, 
which might be interpreted as a disk seen edge-on. The latter com- 
ponent coincides with an edge-on disk observed in maser emission 
jGreenhill et al. I l2003l for the case of the Circinus galaxy). Some 
of these interpretations qualitativ ely match well with t he results of 
our simulations (see discussion in lTristram et al .120071) . It has to be 
taken into account that our modelling is meant to be suitable for an 
average Seyfert galaxy. 



6.1 Observing our standard model 

The basic procedure of our effort of simulating the observational 
properties of dust tori is subdivided into several steps: First of 
all, hydrodynamic simulations are carried out with the help of 
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Figure 8. Images of our standard model at 12 [am. Shown are the indination angles 8°, 30°, 60° 
logarithmic colour scale. Labels are given in pc. 



and 90° for a common azimuthal angle of = 45° with a 



the TRAMP code as described above. In a postprocessing step, 
the resulting gas density distribution is converted into the corre- 
sponding dust density distribution with the help of a constant gas- 
to-dust-ratio of 250 applied to those patches of the gas distribu- 
tion, where temperatures are below the sublimation temperatures 
(Tsub) of the various grain species. Our dust model takes graphite 



1500 K) as well as silicate grains (T|^'b ~ 1000 K) 



,yigraph 
^ sub 

into account. This mixture leads to a prominent spectral feature 
at 9.7 )i.m and a less pronounced feature at 18.5 (xm. A slightly 
higher gas-to-dust rat io compared to the local galactic value of 160 
dSodroski et"ai]| 19941) was applied, as we expect it to increase to- 
wards the centres of active galaxies, as dust is likely to be destroyed 
in the harsh environment close to the energy source. The resulting 
dust density distribut ions are then fed into the 3D radiative transfer 
code MC3D jWolf e t al. 1999; Wolf 2003). Here, in a two-step ap- 
proach, first the temperature distribution is calculated, that is then 
used to simulate observable quantities like spectral energy distribu- 
tions (SEDs) or images at various wavelengths. Grain- size depen- 
dent sublimation (as described in lSchartmann et al]|2008h is only 
considered during the radiative transfer calculations and, therefore, 
only for dust in the close vicinity of the AGN. Here, we split the 
dust opacity model into 6 different grains, made up of the smallest 
(0.0 05 [im) and largest ( 0.25 ^.m) grain oftheMRr^ll- size distribu- 
tion l lMathis etaf]| l977h and the three grain species: astronomical 
silicate and the two orientations of graphite grains (D raine & l3 
ll984lLaor & Drainell993l : IWeingartner & Drainel200ll) . As radia- 
tive transfer calculations are very time and memory consuming, the 
hydrodynamic grid - an assembly of all three domains - has to be 
mapped on to a slightly smaller grid. A resolution of 1 14 grid cells 
in radial direction (instead of 181 for the composed hydrodynamic 
model), 61 grid cells in 9 direction (instead of 94 in the hydro- 
dynamic models) and 128 in t/i-direction (instead of 49 per quad- 
rant in our hydrodynamic models) is currently feasible. We use 96 
distinct wavelengths which are logarithmically distributed between 
0.001 and 2000 [im, with an enhanced resolution in the infrared 
regime and especially around the silicate features. The gas tem- 
perature from the hydrodynamic simulation is only used to decide 



whether dust is present at a certain grid point (Tgas < T. 



sil /graph 
sub 



or not (Tgas ^ T^ub^^ aph^ ^^j. jTg^Q^ jjjg jjygj temper- 

ature during the hydrodynamical simulation explicitly, we set the 



^ The dust model is labelled according to the initials of Mathis, Rumpl and 
Nordsieck. 



dust temperature to zero prior to the radiative transfer calculation 
and assume that the dust distribution is solely heated by the central 
AGN. We assume a bolometric luminosity of the accretion disk of 
Ldisk = 1-2 • 10" Lo, with a cos 6 radiation characteristic, pref- 
erentially emitting perpendicular to the plane of the disk. 

Fig- m shows an inclination angle study for images at 12 [xm 
of the final stage (see Fig.O of our standard model. Shown is the 
dust re-emission and the central radiation of the almost face-on case 
(i = 8°fl i = 30°, i = 60° and the edge-on case (i = 90°) for an 
azimuthal angle of (ji = 45°. Even at an inclination angle of 60°, 
the central source still possesses the highest surface brightness at 
these (ji-angles, but it vanishes behind the dense disk component 
at larger inclination angles (90°). The dense dust disk is the sec- 
ond brightest characteristic for our models. As expected, it does 
not appear smooth, but a filamentary structure is visible, as well as 
a slight nuclear spiral structure (only visible in the high-resolution 
electronic version), caused by the differential rotation of the disk. 
Due to the large optical depth at intermediate inclination angles, 
the torus shows an asymmetric structure with respect to the mid- 
plane, as only the upper, directly illuminated funnel wall is visible. 
The same holds true for the 60° inclination. Here, as well as in the 
edge-on case, shadows of dense - in most cases radially inwards 
moving - filaments are visible. The anisotropic radiation source 
(preferentially emitting along the axis) also contributes to the ver- 
tically elongated bright regions and enhances the dark lanes in the 
midplane. 

In Fig. 121 SEDs of the standard model are presented in terms 
of pure dust re-emission SEDs. Due to the turbulent structure of the 
dense disk component close to the midplane and the filaments and 
blobs of high density gas within or close to the funnel region (com- 
pare to Fig.|6]l, the model for the reduction of the s ilicate feature in 
emission as described in lSchartmann et al] bOOSi) works fine and 
only a moderate emission band with a strength of 0.31 is visible at 
9.7 ^m. For the following analysis, we defi ne the silicate feature 
strength in accordance with lShi et alj ( |2006|) as 



-Ffoat — -Fee 



(6) 



* An inclination angle of 8° is used in order to investigate only lines of 
sight within the computational domain of the hydrodynamical simulations 
as well as within the torus of our previously investigated 2D radiative trans- 
fer simulations we compare to in Sect. 16. 31 
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NGC 4151 (scaled) 
Mrk 841 (scaled) 



100 



Figure 9. SEDs of our standard model. Shown are the inclination angles 
8°, 30°, 60° and 90° (top to bottom) for a common azimuthal angle of 
<A = 45°. 



where Jfoat is the flux at the wavelength of the sihcate feature and 
Jcont is the spline fitted continuum, anchored at wavelengths be- 
tween 5.0 and 7.5 \xm and from 25.0 to 40.0 \im. With increasing 
inclination angle, it changes to moderate absorption: -0.10 for 30°, 
-0.41 for 60° and -0.54 for 90° . This is remarkable, as the same av- 
erage optical depth within the midplane (< rg.r >oqu~ 400!!) 
would cause much deeper absorption - which has never been ob- 
served - for the case of a torus with continuous dust distribution. 
These small feature depths are caused by the fact that the silicate 
emission feature is produced within an elongated volume along the 
funnel region. Thus, the appearance of the silicate feature is due 
to a mixture of emission and absorption and, accordingly, existing 
channels of low absorption towards the observer (though leading to 
silicate emission features) can also contribute to the final morphol- 
ogy of the silicate feature. 

For a detailed discussion of the accuracy and limitations of 
th ese kind of radiative tran sfer calculations, we refer the reader 
to ISchartmann et al.l | |2008|) . The most severe problems occur in 
the midplane, where the largest optical depths are reached, which 
lead to an overestimation of the temperature due to steep gradients. 
Here, a higher spatial resolution would be desirable. However, we 
are confident that our generic results are robust, since a resolution 
study for th e case of SEDs yie l ded al most identical results, as also 
described in lSchartmann et al] ( 1200 Sl) . 



6.2 Comparison with Spitzer spectra 

Figure [TO] shows a comparison of our standard hydrodynamic 
model with observations of two Seyfert galaxies with the Infra-Red 
Spectrograph (IRS) onboard the Spitzer Space Telescope, in high- 
sensitivity, low resolution mode. The solid curve refers to an incli- 
nation angle of 30° towards our standard model. Shown in colour 
are the mid-infrared SEDs of the two intermediate type (Sy 1.5 ac- 
cording to NEeH the NASA/IPAC Ext ragalactic Database) Se yfert 
galaxies NGC 4151 (blue sohd line, IWeedman et all |200^ and 
Mrk 841 (green, H. W. W. Spoon, private communication). Both are 
scaled, in order to obtain fluxes similar to our standard model, as we 
are only interested in a comparison of the shape of the curves. The 
30° -inclination angle of our simulated SED shows a self-absorbed 
silicate feature, yielding a good approximation of the green curve 
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Figure 10. Comparison of our standar d model (solid line - i = 30° ) with 
IRS Spitzer observations of NGC 4151 iWeedman et al.l2005l) and Mrk 841 
(H. W. W. Spoon, private communication). 



in the respective wavelength range. As we are applying pure dust 
radiation transfer, we are not able to produce emission lines in our 
modelled SEDs and we can only compare the underlying contin- 
uum flux. Deviations between observations and model are visible in 
a shift of the maximum of emission of our model towards slightly 
longer wavelengths and a feature around 34 yuca. The latter arises 
most likely from crystalline forsterite in emission (H. W. W. Spoon, 
private communication), which is not included in our dust model. 

As our SEDs do not result from a fitting procedure, this is a 
very noteworthy result. Especially because the physical parameters 
we use are meant to represent a typical Seyfert galaxy. The sharp 
decrease of the flux towards longer wavelengths results from miss- 
ing cold dust in the outer part. 



6.3 Silicate feature strength H i column density relation 

Compared to the previously discussed infrared wavelength range, 
the opacity drops steeply in the X-ray range. Therefore, a certain 
amount of X-ray photons can escape and provides thus the possi- 
bility to determine the hydrogen column density on the line of sight 
observationally. 

For the case of our hydrodynamic models, H I column densi- 
ties are directly integrated along the line of sight. Only cells are 
considered, where the temperature is lower than the ionisation tem- 
perature of hydrogen. However, in our continuous models to which 
we will compare our new hydrodynamical models in the following, 
we cannot take this into account. Therefore, we assume a constant 
gas-to-dust ratio of 250 (see discussion in Sect. 16. Il l, in order to 
calculate column densities from the continuous dust distributions. 
The strength of the silicate feature is derived as given in formula 
[6]and is taken at the wavelength, where the maximum distance be- 
tween the SED and the fitted continuum is reached, as it can be ap- 
preciably shifted for the case of self-absorbed emission. With the 
help of this procedure, our simulati ons can be direct ly co mpared to 
the ob servational data published in lShi et al] ( |2006|) . The lShi et al.l 
1 I2OO6I) sample comprises of 85 AG N of various types, for which 
Spitzer InfraRed Spectrograph (IRS, iHouck et alj 2004h observa- 
tions or ground-based measurements from the lRoche et al. Il ll99ll) 
sample (four objects) exist, as well as X-ray data. H I column den- 
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Figure 11. Comparison of resulting silicate feature strengths and col- 
umn densities of our hydromodels (standard model plus models from a 
supernova and mass loss rate study) given as stars and the continuous 
TTM-models (standard model and models from a dust mass study, see 
ISchartmann et alj|2005h given as triangles. The filled stars denote our hy- 
drodynamical standard model. Colour denotes inclination angle (blue - 8°, 
green - 30° , magenta - 60° , red - 90° , black - intermediate angles of the 
continuous standard model). The so lid lines are linear fits to the observa- 
tional data shown in lShi et aDl200d (black - all objects, magenta - Seyfert 
galaxy sample). The two outlying objects NGC7213 and PG 1351+640 are 
given additionally. Ellipses are drawn to guide the eye. 



sities are obtained from X-ray spectra taken from the Chandra data 
archive directly for 8 objects. For further 77 objects, column den- 
sities were collected from the literature. A correlation was found 
between the silicate feature strength (see equation|6} and the X-ray 
absorption column density in a sense that low H I columns corre- 
spond to silicate emission and high columns correspond to silicate 
absorption: 

Afcaturo = 3.3 ± 0.5 - (0.15 ± 0.02) ■ log(AfHi). (7) 

This is shown as the black line in Fig. [TT] When lstii et al.l j2006l) 
determine the correlation between the silicate feature strength and 
the X-ray colunm density data separately for the various types of 
AGN, a large scatter between them is found, as is also visible in 
the data range for each type of AGN. Indicated in magenta is the 
correlation for Seyfert galaxies alone, given by: 

^feature — 2.6 ± 0.7 - (0.12 ± 0.03) ■ log(ArHi). (8) 

In Fig. [TT] these correlations are compared to our simula- 
tions. Colours of the symbols indicate inclination angle: blue - 
8°, green - 30°, magenta - 60° and red - 90°. For the case of 
inclination angles without any dust along the line-of-sight, we as- 
signed a value of A^^hi = 10^" cm~^, taking foreground gas into 
account. Our hydrodynamic models are given as stars and the tri- 
angles denote results fro m our continuous TTM-model simulations 
jSchartmann et alj2005l) . Overlayed ellipses roughly mark the dis- 
tribution of symbols for the two classes of simulations. Excluding 



well with the observed correlations, whereas the continuous mod- 
els show a much too steep dependence. Although the emission 
features in the Seyfert 1 case are rather too pronounced, they are 
in good agreement with the two outlying galaxies NGC7213 and 
PG 1 35 1-1-640. Furthermore, the silicate feature changes from emis- 
sion to absorption within a very small range of column densi- 
ties and it is also evident that column densities much higher than 
10^'* cm"^ cannot be reached, as they would cause unphysically 
hi gh absorption depths of the silicate feature (as already discussed 



Schartmann et al .l2008h . A further concentration of the dust den- 



the values at 10 cm , which are only due to foreground ex- 
tinction, one can see that the hydrodynamical models agree quite 



sity towards the minimum of the potential within the midplane 
might cause a flattening of this distrib ution, but th e n, too narrow 
infrared bumps will be produced, see ISch artmann ('2003'). Obvi- 
ously, the observations presented in lShi et al.. (2006.) are in favour 
of clumpy or filamentary models instead of continuous matter dis- 
tributions, concerning the models we investigated so far. At least, 
models are needed which feature high column densities in combi- 
nation with moderate emission as well as absorption features. In 
the hydrodynamical models presented in this paper, this is possible 
due to the clumpy and filamentary nature of the dust distribution, 
which enables the coexistence of almost free lines of sight towards 
the inner region, where the silicate feature in emission is produced, 
together with highly obscured lines of sight. The highest column 
densities are reached for lines of sight through the dense and tur- 
bulent disk in our models. Furthermore, a remarkable consequence 
of this is the very tight linear correlation of the continuous models 
and the broad distribution for the case of our (clumpy) hydrody- 
namical models. Therefore, the broad scatter in the observational 
d ata might really b e taken as a case for clumpy models, as claimed 
in lShi et al] l l2006l) . The observed s hallow slope of the disc ussed re- 
lation has also been investigated bv lLevenson et alj l l2007h on basis 
of radiative transfer models of simplified geometries, representing 
dust clouds and continuous dust distributions. Whereas the latter 
can yield arbitrarily deep silicate absorption features in their simu- 
lations, single dust clouds can not, as long as their dimensions are 
much smaller than their distance to the illuminating source. With 
the help of these simulations they conclude that the tori of Seyfert 
galaxies must be clumpy, as their mean SED shows only moder- 
ate emission as well as absorption features and ULIRGs (Ultra- 
luminous Infrared Galaxies) must possess nuclear sources, which 
are deeply embedded into a geometrically and optically thick and 
smooth distributi on of dust, as th eir mean SED features deep sil- 
icate absorption l lHao et al.|[2007h . Similar conclusions have been 
drawn with the help of a c lumpy torus model already discussed in 
sectionfTlbv iNenkova et alj (I2002i) and .Nenkova et al.. (2008a .h). 



7 SUMMARY 

We present a physical model for the evolution of gaseous tori in 
Seyfert galaxies, in which the stellar evolution of a young nuclear 
star cluster is predominantly responsible for setting up a torus-like 
dust (and gas) distribution. The star cluster provides both, (i) mass 
injection, which is mainly due to planetary nebulae and (ii) injec- 
tion of energy from supernova explosions. Within the course of 
the simulations, a highly dynamical multiphase medium evolves, 
which is comparable to high resol ution models of the supernova- 
driven turbulent ISM in g alaxies tBreitschwerdt & de A villez 200^; 
Ijoung & Mac Lowll2006l) . There, blast waves from supernovae ex- 
plosions sweep up the ISM by transferring kinetic energy to the 
gas, which leads to the subsequent formation of thin shells. These 
shells interact with one another and with already existing filaments. 
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Radiative cooling then produces cold clumps. In our scenario, the 
ISM is constantly refilled by stellar mass loss processes distributed 
according to the stars in the nuclear star cluster. The main structur- 
ing agent of the gas and dust in the torus is the interplay between 
mass input, supernova blast waves locally sweeping away the ISM, 
the streams of cold gas towards the centre as well as cooling insta- 
bility. 

Within this scenario, cold dense clouds and filaments, where dust 
can exist, form naturally. The implemented optically thin cooling 
reduces thermal pressure and, therefore, material streams towards 
the minimum of the potential. This forms long radial filaments, 
with lots of substructure due to interac tion processes. In agreement 
with M HD simulations of the ISM by ISreitschwerdt & de Aville j 
( |2007|) . we find absence of pressure equilibrium of the temperature 
phases. This is the case, as mixing due to supersonic motions and 
turbulence as well as cooling processes act on shorter time-scales 
than relaxation processes. Although no pressure equilibrium can 
be reached, the flow is in a global dynamical equilibrium state. The 
characteristics of this equilibrium delicately depend on the driving 
parameters - that are energy input (supernova rate) and mass in- 
put rate, as well as the gravitational potential, the cooling curve 
and rotation. We find completely different velocities and streaming 
patterns for the various temperature phases. Whereas cold material 
generally tends to move towards the centre, the radial velocity of 
the hot ionised medium shows both inflow (in the inner region) and 
outflow (further out). Depending on the ratio between supernova 
rate and mass loss rate, the inflow/outflow transition shifts towards 
smaller radii for higher SN-to-mass injection ratios and towards 
larger radii, when mass input dominates over supernova heating. 
For very high supernova rates, we observe solely outflow solutions 
for the hot gas component. This solution is likely to dominate the 
initial supernova type II phase shortly after the birth of the stellar 
population. 



The resulting density and temperature structure of the gas is 
subsequently fed into the radiative transfer code MC3D in order to 
obtain observable quantities (SEDs and surface brightness distribu- 
tions in the mid-infrared), which are then compared to recent obser- 
vations. A comparison with spectral energy data from the Spitzer 
satellite shows that, despite of its severe limitations in resolution 
and astrophysics considered within this exploratory study, our hy- 
drodynamical standard model describes one class of Spitzer SEDs 
fairly well, without the necessity of any finetuning of parameters. 
For example the two intermediate type Seyfert galaxies NGC4151 
as well as Mrk 841 belong to this class (see Fig.llOt. The success of 
the comparison with the silicate feature strength to Hi column den- 
sity relation for the hydrodynamical models and the failure of the 
continuous TTM-models leads to the important conclusion that the 
presence of a disk-like density enhancement within the equatorial 
plane in combination with a fluffy structure surrounding it is needed 
in order to get agreement with observational data. Otherwise, large 
X-ray column densities cannot simultaneously be obtained together 
with moderate silicate absorption feature strengths at viewing an- 
gles within the midplane. The latter are guaranteed by a mixture 
of almost dust-free and heavily obscured lines of sight towards the 
hot inner part of the torus close to the funnel. Therefore, we see a 
much weaker silicate absorption feature. Due to these remarkable 
similarities of our exploratory study and observations, we feel en- 
couraged to implement more physical effects into our simulations 
and to study the central gas and dust structure over a longer evolu- 
tionary time of the nuclear star cluster. 
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